# Replication file for:

# The Effect of Wartime Legacies on Electoral Mobilization after Civil War
# Felix Haass, Martin Ottmann

# Datestamp: 29.09.2021

# Description: This file replicates all tables and figures in the main paper

# Usage: please follow the instructions in the file README.txt


# Check if packages are installed -----------------------------------------

list.of.packages <- c("here",
                      "lfe",
                      "fastDummies",
                      "tidyverse",
                      "patchwork",
                      "fixest",
                      "lmtest", 
                      "multiwayvcov",
                      "here",
                      "RcppArmadillo",
                      "Rcpp",
                      "geosphere",
                      "ggsn",
                      "stargazer",
                      "kableExtra",
                      "sf",
                      "data.table",
                      "ggpubr",
                      "modelsummary")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)


# Load Libraries ---------------------------------------------------------------

library(here)
library(lfe)
library(fastDummies)
library(tidyverse)
library(patchwork)
library(fixest)
library(here)
library(lmtest)
library(multiwayvcov)
library(stargazer)
library(kableExtra)
library(sf)
library(ggpubr)
library(ggsn)
library(modelsummary)
library(RcppArmadillo)
library(Rcpp)
library(geosphere)
library(data.table)


# Load data ---------------------------------------------------------------

# RAN data
model_data_kec_month <- readRDS(here("./data/model_data_kec_month.rds"))
model_data_kab_month <- readRDS(here("./data/model_data_kab_month.rds"))

# survey data 
vh_clean <- readRDS(here("./data/vh_clean.rds"))

# geometries for map outlines
aceh_gadm_crosswalk_kec_2005 <- readRDS(here("./data/aceh_gadm_crosswalk_kec_2005.rds"))
aceh_gadm_crosswalk_kab_2005 <- readRDS(here("./data/aceh_gadm_crosswalk_kab_2005.rds"))
sumatera_utara <- readRDS(here("./data/sumatera_utara.rds"))
aceh_gadm_kec <- readRDS(here("./data/aceh_gadm_kec.rds"))


# Figure 1 ----------------------------------------------------------------

## GAM support map ----------------------------------------------------------

gam_control <- ggplot(aceh_gadm_crosswalk_kec_2005) + 
  
  # basic fill 
  geom_sf(aes(fill = control_nvms), 
          color = "darkgrey") +
  
  # add kabupaten outline
  geom_sf(data= aceh_gadm_crosswalk_kab_2005, alpha = 0, 
          color = "black") + 
  scale_fill_distiller("GAM support", palette = "YlOrRd", direction = 1) +
  
  # scalebar + north
  north(aceh_gadm_crosswalk_kec_2005 %>% select(control_nvms)) + 
  scalebar(data = aceh_gadm_crosswalk_kec_2005 %>% select(control_nvms), 
           dist = 50, dist_unit = "km",
           transform = TRUE, model = "WGS84", location = "bottomleft", st.size = 4) +
  
  # add Northern Sumatra outline
  geom_sf(data = sumatera_utara %>%
            st_transform(4326) %>%
            st_crop(st_bbox(aceh_gadm_crosswalk_kec_2005)),
          color = "darkgrey", fill = "lightgrey") +
  
  # # outline of Aceh
  geom_sf(data = aceh_gadm_kec  %>% st_union(),
          color = "black",  alpha = 0) +
  theme_bw() +
  coord_sf(datum = NA) +
  labs(x = "", y = "", title = "") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        # panel.border = element_blank(),
        axis.text = element_blank(), 
        legend.position = "bottom") 


## Tsunami aid map ---------------------------------------------------------

tsunami_aid_map_data <- model_data_kec_month %>% 
  mutate(ID2005 = as.character(ID2005)) %>% 
  group_by(ID2005) %>% 
  summarise(number_of_new_projects_all = sum(number_of_new_projects, na.rm = T), 
            number_of_new_projects_pre2006 = sum(number_of_new_projects[year <= 2006], na.rm = T),
            number_of_new_projects_post2006 = sum(number_of_new_projects[year > 2006], na.rm = T),) %>% 
  left_join(., aceh_gadm_crosswalk_kec_2005,
            by = "ID2005") %>% 
  st_as_sf()

tsunami_aid_map <- ggplot(tsunami_aid_map_data) + 
  # basic fill 
  geom_sf(aes(fill = log(number_of_new_projects_all + 1)), 
          color = "darkgrey") +
  # add kabupaten outline
  geom_sf(data= aceh_gadm_crosswalk_kab_2005,  alpha = 0, 
          color = "black") + 
  scale_fill_distiller("Number of projects (log + 1)", palette = "PuBu", na.value = "green", 
                       direction = 1) +
  # scalebar + north
  north(aceh_gadm_crosswalk_kec_2005 %>% select(control_nvms)) + 
  scalebar(data = aceh_gadm_crosswalk_kec_2005 %>% select(control_nvms), 
           dist = 50, dist_unit = "km",
           transform = TRUE, model = "WGS84", location = "bottomleft", st.size = 4) +
  # add Northern Sumatra outline
  geom_sf(data = sumatera_utara %>%
            st_transform(4326) %>%
            st_crop(st_bbox(aceh_gadm_crosswalk_kec_2005)),
          color = "darkgrey", fill = "lightgrey") +
  # # outline of Aceh
  geom_sf(data = aceh_gadm_kec  %>% st_union(),
          color = "black",  alpha = 0) +
  theme_bw() +
  coord_sf(datum = NA) +
  labs(x = "", y = "", title = "") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        # panel.border = element_blank(),
        axis.text = element_blank(), 
        legend.position = "bottom") 


## Scatterplot -------------------------------------------------------------

aid_control_scatterplot <- ggplot(tsunami_aid_map_data, 
                                  aes(x = log(control_nvms +1), 
                                      y= log(number_of_new_projects_all +1))) + 
  geom_point(size = 3, alpha = 0.5) + 
  geom_smooth(method = "lm", se = F, size = 2) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        # panel.border = element_blank(),
        # axis.text = element_blank(), 
        legend.position = "bottom") +
  labs(x = "GAM support (log + 1)", y = "Number of new projects (log + 1)")


## Output final figure 1 ---------------------------------------------------

ggsave(gam_control + tsunami_aid_map + aid_control_scatterplot + plot_annotation(tag_levels = "A"), 
       height = 7, 
       width = 14, scale = 0.8, filename = here("./figures/paper_figure_1.pdf"))


# Table 1  -------------------------------------------------------------

## Panel A: RAN data models -------------------------------------------------------

# kecamatan
aid_model_kec_fe_time <- felm(log(number_of_new_projects + 1) ~ 
                                pre_election_2006:control_dummy_median_nvms 
                              |  kecid + yearmonth_fct | 0 | kecid , 
                              data = model_data_kec_month )

# kabubaten
aid_model_kab_fe_time <- felm(log(number_of_new_projects + 1) ~ 
                                pre_election_2006:control_dummy_median_nvms
                              |  KABNO + yearmonth | 0 | KABNO , 
                              data = model_data_kab_month)

# USD as DV
# kecamatan
aid_model_kec_fe_time_usd <- felm(log(committed_usd + 1) ~
                                    pre_election_2006:control_dummy_median_nvms
                                  
                                  |  kecid +yearmonth_fct | 0 | kecid,
                                  data = model_data_kec_month )

# kabupaten
aid_model_kab_fe_time_usd <- felm(log(committed_usd + 1) ~ 
                                    pre_election_2006:control_dummy_median_nvms
                                  |  KABNO + yearmonth | 0 | KABNO , 
                                  data = model_data_kab_month)

# LPM
aid_model_lpm <- felm(I(number_of_new_projects > 0) ~ 
                        pre_election_2006:control_dummy_median_nvms
                      |  kecid + yearmonth_fct | 0 | kecid , 
                      data = model_data_kec_month  )


### Output RAN data models table --------------------------------------------

stargazer::stargazer(aid_model_kec_fe_time, 
                     aid_model_kab_fe_time,
                     aid_model_kec_fe_time_usd,
                     aid_model_kab_fe_time_usd,
                     aid_model_lpm,
                     out = here("./tables/paper_table_1_panel_A.html"),
                     covariate.labels = c("GAM Area X Pre-Election") ,
                     dep.var.caption = "", 
                     dep.var.labels = c("No. of New Projects (log +1)",
                                        "Aid Volume (log + 1)",
                                        "Pr(New Aid Project > 0)"),
                     # omit = c(3,4,6),
                     omit = "as.numeric|Constant",
                     add.lines = list( c("Geographical Unit", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}"),
                                       c("Unit FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                       c("Month FE", rep("\\multicolumn{1}{c}{Yes}", 8))),
                     float = F, 
                     digits = 2,
                     notes = "",
                     notes.append = F, 
                     notes.label = "", 
                     align = T, 
                     keep.stat = c("n", "adj.rsq"))


## Panel B: Survey data models ------------------------------------------------------

### Non-election year 2005 --------------------------------------------

baseline_lpm_05 <- felm(election_tsunami_05 ~
                          gam_support 
                        | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_lpm_05 <- felm(election_tsunami_05 ~
                           gam_support +
                           # structure
                           H060i_grp_kpa  +
                           # demographics
                           I(tsunami_affected > 0) +
                           log(H027_vill_pop) +
                           log(H032a_perc_aceh+1) +
                           factor(H036_idp) +
                           H048_poorHH +
                           H053a_primary98 +
                           number_of_vh +
                           # factor(H053f_mos98) +
                           I(H110_conf_victims0105 > 0) +
                           # interviewer FE
                           factor(H147_lang) +
                           factor(Hq152_r_6_KPA)
                         | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)


baseline_logit_05 <- glm(election_tsunami_05 ~
                           gam_support + 
                           factor(H001a_regency),
                         family="binomial", 
                         data = vh_clean)

baseline_logit_se_05 <- lmtest::coeftest(baseline_logit_05, 
                                         multiwayvcov::cluster.vcov(baseline_logit_05, 
                                                                    vh_clean$H001c_subdistrict))
baseline_logit_se_05 <- broom::tidy(baseline_logit_se_05)$std.error



covariate_logit_05 <- glm(election_tsunami_05 ~
                            gam_support + 
                            H060i_grp_kpa  +
                            # demographics
                            I(tsunami_affected > 0) +
                            log(H027_vill_pop) +
                            log(H032a_perc_aceh+1) +
                            factor(H036_idp) +
                            H048_poorHH +
                            H053a_primary98 +
                            number_of_vh +
                            # factor(H053f_mos98) +
                            I(H110_conf_victims0105 > 0) +
                            # interviewer FE
                            factor(H147_lang) +
                            factor(Hq152_r_6_KPA) +
                            factor(H001a_regency),
                          family="binomial", 
                          data = vh_clean)

covariate_logit_se_05 <- lmtest::coeftest(covariate_logit_05, 
                                          multiwayvcov::cluster.vcov(covariate_logit_05, 
                                                                     vh_clean$H001c_subdistrict))

covariate_logit_se_05 <- broom::tidy(covariate_logit_se_05)$std.error


### Election_year 2006 -------------------------------------------------

baseline_lpm <- felm(election_tsunami_06 ~
                       gam_support 
                     | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_lpm <- felm(election_tsunami_06 ~
                        gam_support +
                        # structure
                        H060i_grp_kpa  +
                        # demographics
                        I(tsunami_affected > 0) +
                        log(H027_vill_pop) +
                        log(H032a_perc_aceh+1) +
                        factor(H036_idp) +
                        H048_poorHH +
                        H053a_primary98 +
                        number_of_vh +
                        I(H110_conf_victims0105 > 0) +
                        # interviewer FE
                        factor(H147_lang) +
                        factor(Hq152_r_6_KPA)
                      | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

baseline_logit <- glm(election_tsunami_06 ~
                        gam_support + 
                        factor(H001a_regency),
                      family="binomial", 
                      data = vh_clean)

baseline_logit_se <- lmtest::coeftest(baseline_logit, 
                                      multiwayvcov::cluster.vcov(baseline_logit, 
                                                                 vh_clean$H001c_subdistrict))
baseline_logit_se <- broom::tidy(baseline_logit_se)$std.error

covariate_logit <- glm(election_tsunami_06 ~
                         gam_support + 
                         H060i_grp_kpa  +
                         # demographics
                         I(tsunami_affected > 0) +
                         log(H027_vill_pop) +
                         log(H032a_perc_aceh+1) +
                         factor(H036_idp) +
                         H048_poorHH +
                         H053a_primary98 +
                         number_of_vh +
                         # factor(H053f_mos98) +
                         I(H110_conf_victims0105 > 0) +
                         # interviewer FE
                         factor(H147_lang) +
                         factor(Hq152_r_6_KPA) +
                         factor(H001a_regency),
                       family="binomial", 
                       data = vh_clean)

covariate_logit_se <- lmtest::coeftest(covariate_logit, 
                                       multiwayvcov::cluster.vcov(covariate_logit, 
                                                                  vh_clean$H001c_subdistrict))

covariate_logit_se <- broom::tidy(covariate_logit_se)$std.error


### Output survey models table  ------------------------------------------------------

stargazer::stargazer(baseline_lpm_05, 
                     covariate_lpm_05,
                     baseline_logit_05, 
                     covariate_logit_05, 
                     baseline_lpm, 
                     covariate_lpm,  
                     baseline_logit, 
                     covariate_logit,
                     out = here("./tables/paper_table_1_panel_B.html"),
                     type = "latex", 
                     dep.var.caption = "",
                     dep.var.labels = c("Did village receive tsunami aid in 2005?", 
                                        "Did village receive tsunami aid in 2006?"),
                     covariate.labels = c("Majority in village supported GAM"),
                     omit = "H0|vh|tsunami|lang|Hq|Constant|H110",
                     add.lines = list(c("Kab. FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                      c("Covariates", rep(c("\\multicolumn{1}{c}{No}", 
                                                            "\\multicolumn{1}{c}{Yes}"), 4))), 
                     # notes = strwrap(c("Linear probability models. Robust standard errors, clustered by Kecamatan in parentheses. Covariates: Tsunami affected, Village Pop. (Q29), Percent Ethnic Acehnese (Q32), IDPs (Q36), Poverty (Q48), Primary Schools (Q53), Conflict Victims (Q110), Interviewer Language - Acehnese (Q147),Interviewer Language - Other (Q147),KPA present during interview (Q152)"), 70),  
                     float = F, 
                     digits = 2,
                     notes = "",
                     se = list(NULL, NULL, baseline_logit_se, covariate_logit_se),
                     notes.append = F, 
                     notes.label = "", 
                     align = T, 
                     model.names = F,
                     column.labels = c("OLS", "Logit", "OLS", "Logit"),
                     column.separate = c(2, 2, 2, 2), 
                     keep.stat = c("n", "adj.rsq", "aic"))


# Figure 2: Organizational legacies ----------------------------------------

source(here("./code/functions/interaction_plots_felm.R"))

baseline_kpa_lpm <- felm(election_tsunami_06 ~
                           gam_support * 
                           H060i_grp_kpa  
                         | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)


covariate_kpa_lpm <- felm(election_tsunami_06 ~
                            gam_support *
                            H060i_grp_kpa +
                            # demographics
                            I(tsunami_affected > 0) +
                            H029_vill_pop98 +
                            log(H032a_perc_aceh+1) +
                            factor(H036_idp) +
                            H048_poorHH +
                            H053a_primary98 +
                            I(H110_conf_victims0105 > 0) +
                            # interviewer FE
                            factor(H147_lang) +
                            factor(Hq152_r_6_KPA)
                          | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)


covariate_youth_lpm <- felm(election_tsunami_06 ~
                              gam_support *
                              H060g_grp_youth +
                              # demographics
                              I(tsunami_affected > 0) +
                              H029_vill_pop98 +
                              log(H032a_perc_aceh+1) +
                              factor(H036_idp) +
                              H048_poorHH +
                              H053a_primary98 +
                              I(H110_conf_victims0105 > 0) +
                              # interviewer FE
                              factor(H147_lang) +
                              factor(Hq152_r_6_KPA)
                            | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_ethnic_lpm <- felm(election_tsunami_06 ~
                               gam_support *
                               H060e_grp_cult +
                               # demographics
                               I(tsunami_affected > 0) +
                               H029_vill_pop98 +
                               log(H032a_perc_aceh+1) +
                               factor(H036_idp) +
                               H048_poorHH +
                               H053a_primary98 +
                               I(H110_conf_victims0105 > 0) +
                               # interviewer FE
                               factor(H147_lang) +
                               factor(Hq152_r_6_KPA)
                             | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_women_lpm <- felm(election_tsunami_06 ~
                              gam_support *
                              H060h_grp_women +
                              # demographics
                              I(tsunami_affected > 0) +
                              H029_vill_pop98 +
                              log(H032a_perc_aceh+1) +
                              factor(H036_idp) +
                              H048_poorHH +
                              H053a_primary98 +
                              I(H110_conf_victims0105 > 0) +
                              # interviewer FE
                              factor(H147_lang) +
                              factor(Hq152_r_6_KPA)
                            
                            | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_pol_lpm <- felm(election_tsunami_06 ~
                            gam_support *
                            H060f_grp_pol +
                            # demographics
                            I(tsunami_affected > 0) +
                            H029_vill_pop98 +
                            log(H032a_perc_aceh+1) +
                            factor(H036_idp) +
                            H048_poorHH +
                            H053a_primary98 +
                            I(H110_conf_victims0105 > 0) +
                            # interviewer FE
                            factor(H147_lang) +
                            factor(Hq152_r_6_KPA)
                          
                          | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

covariate_rel_lpm <- felm(election_tsunami_06 ~
                            gam_support *
                            H060d_grp_rel +
                            # demographics
                            I(tsunami_affected > 0) +
                            H029_vill_pop98 +
                            log(H032a_perc_aceh+1) +
                            factor(H036_idp) +
                            H048_poorHH +
                            H053a_primary98 +
                            I(H110_conf_victims0105 > 0) +
                            # interviewer FE
                            factor(H147_lang) +
                            factor(Hq152_r_6_KPA)
                          
                          | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)


# Plotting marginal effects 

covariate_kpa_me <- interaction_plot_binary(covariate_kpa_lpm, 
                                            "gam_support", 
                                            "H060i_grp_kpa", 
                                            "gam_support:H060i_grp_kpa") %>% 
  mutate(assoc = "KPA", 
         present = c("No", "Yes"))


covariate_ethnic_me <- interaction_plot_binary(covariate_ethnic_lpm, 
                                               "gam_support", 
                                               "H060e_grp_cult", 
                                               "gam_support:H060e_grp_cult") %>% 
  mutate(assoc = "Ethnic", 
         present = c("No", "Yes"))


covariate_women_me <- interaction_plot_binary(covariate_women_lpm, 
                                              "gam_support", 
                                              "H060h_grp_women", 
                                              "gam_support:H060h_grp_women") %>% 
  mutate(assoc = "Women", 
         present = c("No", "Yes"))


covariate_pol_me <- interaction_plot_binary(covariate_pol_lpm, 
                                            "gam_support", 
                                            "H060f_grp_pol", 
                                            "gam_support:H060f_grp_pol") %>% 
  mutate(assoc = "Political", 
         present = c("No", "Yes"))


covariate_rel_me <- interaction_plot_binary(covariate_rel_lpm, 
                                            "gam_support", 
                                            "H060d_grp_rel", 
                                            "gam_support:H060d_grp_rel") %>% 
  mutate(assoc = "Religious", 
         present = c("No", "Yes"))


covariate_youth_me <- interaction_plot_binary(covariate_youth_lpm, 
                                              "gam_support", 
                                              "H060g_grp_youth", 
                                              "gam_support:H060g_grp_youth") %>% 
  mutate(assoc = "Youth", 
         present = c("No", "Yes"))


all_marginal_effects <- bind_rows(covariate_kpa_me, 
                                  covariate_youth_me, 
                                  covariate_pol_me, 
                                  covariate_rel_me, 
                                  covariate_women_me, 
                                  covariate_ethnic_me) %>% 
  group_by(assoc) %>% 
  mutate(positive_effect = beta[present == "Yes"]) %>% 
  ungroup() %>% 
  mutate(present = forcats::fct_relevel(present, c("Yes", "No"))) %>% 
  mutate(assoc = forcats::fct_reorder(assoc, positive_effect * -1)) # order labels in reverse order


# build ME plot

marginal_effects_assoc <- ggplot(all_marginal_effects, aes(x = present, y = beta)) + 
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = se_low, ymax = se_up), 
                width = 0) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  facet_wrap(~ assoc ) +
  labs(x = "\nIs an association of type ... present in village?", 
       y = "\nProbabilty of Receiving\nTsunami Aid in 2006") +
  theme_bw() + 
  theme(panel.grid = element_blank()) 

ggsave(marginal_effects_assoc, filename = here("./figures/paper_figure_2.pdf"), 
       width = 9, height = 6, scale = 0.8)


# Figure 3: Oversight & Discretionary control ----------------------------------------------------------------

partner_private_count_model <- felm(log(partner_private_count + 1) ~ 
                                      pre_election_2006:control_dummy_median_nvms 
                                    |  kecid +yearmonth_fct | 0 | kecid , 
                                    data = model_data_kec_month )


partner_gov_count_model <- felm(log(partner_gov_count + 1) ~ 
                                  pre_election_2006:control_dummy_median_nvms
                                |  kecid +yearmonth_fct | 0 | kecid , 
                                data = model_data_kec_month )


partner_io_count_model <- felm(log(partner_io_count + 1) ~ 
                                 pre_election_2006:control_dummy_median_nvms
                               |  kecid +yearmonth_fct | 0 | kecid , 
                               data = model_data_kec_month )

partner_natl_ngo_count_model <- felm(log(partner_natl_ngo_count + 1) ~ 
                                       pre_election_2006:control_dummy_median_nvms
                                     |  kecid +yearmonth_fct | 0 | kecid , 
                                     data = model_data_kec_month )

partner_intl_ngo_count_model <- felm(log(partner_intl_ngo_count + 1) ~ 
                                       pre_election_2006:control_dummy_median_nvms
                                     |  kecid +yearmonth_fct | 0 | kecid , 
                                     data = model_data_kec_month )

partner_bilateral_count_model <-  felm(log(partner_bilateral_count + 1) ~ 
                                         pre_election_2006:control_dummy_median_nvms
                                       |  kecid +yearmonth_fct | 0 | kecid , 
                                       data = model_data_kec_month )


partner_international_private_count_model <-  felm(log(partner_international_private_count + 1) ~ 
                                                     pre_election_2006:control_dummy_median_nvms
                                                   |  kecid +yearmonth_fct | 0 | kecid , 
                                                   data = model_data_kec_month )


partner_models <- bind_rows(broom::tidy(partner_private_count_model) %>% mutate(model = "Private Sector"),
                            broom::tidy(partner_gov_count_model) %>% mutate(model = "Government Agency"),
                            broom::tidy(partner_international_private_count_model) %>% mutate(model = "International Private"),
                            broom::tidy(partner_io_count_model) %>% mutate(model = "UN"),
                            broom::tidy(partner_natl_ngo_count_model) %>% mutate(model = "National NGO"),
                            broom::tidy(partner_intl_ngo_count_model) %>% mutate(model = "International NGO"),
                            broom::tidy(partner_bilateral_count_model) %>% mutate(model = "Bilateral"))

partner_models <- partner_models %>% 
  mutate(model = str_wrap(model, 10)) %>% 
  mutate(model = forcats::fct_reorder(model, estimate * -1))

coefplot_partner <- ggplot(partner_models, aes(x = model, y = estimate)) +
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error, 
                    ymax = estimate + 1.96 * std.error), 
                width = 0) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Partner Agency", y = "Coefficient Estimate") +
  # coord_flip() +
  theme_bw() + 
  theme(panel.grid = element_blank()) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

ggsave(coefplot_partner, filename = here("./figures/paper_figure_3.pdf"), 
       height = 4, 
       width = 11, 
       scale = 0.8)
